function Moments = ComputeMoments( Params , Grids, L , Lh , zeta , v , GE )

run MiscFunctions.m

% Fraction of population
w = L .* Grids.mx ; w = w / sum(w) ;
nw = Grids.mx ; nw = nw / sum(nw) ;
w0 = ones(size(zeta)); w0 = w0 / sum(w0) ;

% Unemployment aggregates
ux = uRate(v,zeta) ;

    % Weighted
MeanU = sum( ux .* w ) ;
StdU = sqrt( sum( ux.^2 .* w ) - MeanU^2 ) ;
    % Unweighted
MeanU0 = sum( ux .* w0 ) ;
StdU0 = sqrt( sum( ux.^2 .* w0 ) - MeanU0^2 ) ;


% Variance decomposition: weighted
lsep  = log( Params.mu + Params.delta * zeta ) ;
lfind = log( Phi(v,zeta) ) ;

vSep       = Variance( lsep , w ) ;
vFind      = Variance( lfind , w ) ;
covSepFind = - 2 * Covariance( lsep , lfind , w ) ;
vU         = Variance( lsep - lfind , w ) ;
vUfull     = Variance( log( ux ./ ( 1 - ux ) ) , w ) ;

% Variance decomposition: unweighted
vSep0       = Variance( lsep , w0 ) ;
vFind0      = Variance( lfind , w0 ) ;
covSepFind0 = - 2 * Covariance( lsep , lfind , w0 ) ;
vU0         = Variance( lsep - lfind , w0 ) ;
vU0full     = Variance( log( ux ./ ( 1 - ux ) ) , w0 ) ;

% Value of moving M0
U2x =  Grids.x.^( Params.eps * ( 1 + Params.eta - Params.omega ) / ( Params.omega + Params.eps * ( 1 + Params.eta ) ) ) ...
   .* ( Params.b + v ).^( Params.eps * ( 1 + Params.eta + Params.psi ) / Params.P1 ) ...
   .* G(v,zeta).^( - Params.eps * ( Params.omega + Params.psi ) / Params.P1 ) ;

px = (   1 / ( Params.omega + Params.eps * ( 1 + Params.eta ) )  ...
       * ( 1 + Params.eta * Params.psi / Params.P1 ) )  ...
      * Params.aCondx ;
  
p2x = (  Params.eps / ( Params.omega + Params.eps * ( 1 + Params.eta ) )  ...
       * ( 1 + Params.eta * Params.psi / Params.P1 ) )  ...
      * Params.aCondx ;
  
Dtildex =   exp( Params.stdaCondx^2 / 2 * px^2 ) ...
          * Grids.x.^( px * Params.aCondx ) ;
      
Dx      =   exp( Params.stdaCondx^2 / 2 * p2x^2 ) ...
          * Grids.x.^( p2x * Params.aCondx ) ;
      
% tax rate
tau0 = Tau0(v,zeta) ;
I    = AgLabIncome(v,zeta) ;
E    = AgTaxExp(tau0,Grids.tax,v,zeta,Grids) ;
xii  = E / I ;

% M0: two equivalent calculations (up to numerical approximation)
M0    = ( sum( U2x.^Params.nu .* Dtildex .* Grids.mx ) ).^( Params.P1 / ( 1 + Params.eta + Params.psi ) ) ;
M0alt = (   (1-xii)^(-(Params.omega+Params.psi)/Params.P1) ...
          * sum( Lh .* Grids.mx ) ...
        )^( Params.P1 / ( 1 + Params.eta + Params.psi ) ) ;

% Human capital
x0 = Params.muLearn / ( Params.muLearn + Params.DeltaLearn ) ;                
HC = MeanHC(x0,Lh,v,zeta).^( 1 / Params.P2 ) ;

hx = HC * Params.DLearn ./ ( Params.DLearn + Params.phi * ux ) ;

% aggregate welfare
JM0 = M0^Params.P2 ;
ES = Es(zeta) ;

Welfare = ( 1 - xii) * M0 * sum( ( 1 + Params.barg * ( 1 - ux ) .* ES ) .* hx .* w ) ;

WM = ( 1 - xii) * M0 ;
WE = sum( ( 1 + Params.barg * ( 1 - ux ) .* ES ) .* w ) ;
WH = sum( ( 1 + Params.barg * ( 1 - ux ) .* ES ).* hx .* w ) / WE ;

% local welfare
WUx = ( 1 - xii )^Params.P2 * JM0 * U2x .* Dx ;
WEx = ( 1 + Params.barg * ( 1 - ux ) .* ES ) ;
WHx = hx ;
Wx  = WUx .* WEx .* WHx ;

% entry
j0     = (1-Params.barg)*Params.barg^(Params.alpha/(1-Params.alpha))/Params.rho;
GECnet = GE.C / tau0^Params.gamma ;
Moments.Me     = ...
  GECnet / sum( Lh .* Grids.mx ) ...
         / (   M0^( Params.psi / Params.P1 )^Params.gamma ...
             * MeanHC(x0,Lh,v,zeta)^(1/Params.P2) ...
             * j0 ) ;
       
% aggregate profits and land rents
Moments.TotalProfits   = Moments.Me * AgProfits(v,zeta,Grids) ;

LocalLandRents = (1-xii)^(1-Params.eps*Params.psi/Params.P1) ...
               * M0^(1-Params.psi/Params.P1) ...
               * Grids.x .* LocalHC(v,zeta).^(1/Params.P2) ...
              .* ( Params.b + v ).^(-Params.psi/Params.P1) ...
              .* G(v,zeta).^( 1 - Params.eps * Params.psi / Params.P1 ) ...
              .* Lh / sum(Lh.*Grids.mx) ;

Moments.TotalLandRents = sum( LocalLandRents .* Grids.mx ) ;

% assign output
Moments.vSep        = vSep ;
Moments.vFind       = vFind ;
Moments.covSepFind  = covSepFind ;
Moments.vSep2       = vSep + covSepFind / 2;
Moments.vFind2      = vFind + covSepFind / 2 ;
Moments.vSep2p      = 100 * Moments.vSep2 / vU ;
Moments.vFind2p     = 100 * Moments.vFind2 / vU ;

Moments.vU          = vU ;
Moments.MeanU       = MeanU ;
Moments.StdU        = StdU ;

Moments.vSep0       = vSep0 ;
Moments.vFind0      = vFind0 ;
Moments.covSepFind0 = covSepFind0 ;
Moments.vSep02      = vSep0 + covSepFind0 / 2;
Moments.vFind02     = vFind0 + covSepFind0 / 2 ;
Moments.vSep02p     = 100 * Moments.vSep02 /vU0 ;
Moments.vFind02p    = 100 * Moments.vFind02 /vU0 ;

Moments.vU0         = vU0 ;
Moments.MeanU0      = MeanU0 ;
Moments.StdU0       = StdU0 ;

Moments.Welfare     = Welfare ;
Moments.M0          = M0 ;
Moments.Wx          = Wx ;
Moments.WUx         = WUx ;
Moments.WEx         = WEx ;
Moments.WHx         = WHx ;
Moments.hx          = hx ;
Moments.HC          = HC ;
Moments.U2x         = U2x ;
Moments.ES          = ES ;
Moments.WM          = WM ;  
Moments.WE          = WE ;
Moments.WH          = WH ;
Moments.vUfull      = vUfull ;
Moments.vU0full     = vU0full ;

Moments.M0alt       = M0alt ;

Moments.tau0 = tau0 ;
Moments.I    = I ;
Moments.E    = E ;
Moments.xi   = xii ;


end

